home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / poch.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  13.3 KB  |  436 lines
  1. /* specfunc/poch.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_log.h>
  27. #include <gsl/gsl_sf_pow_int.h>
  28. #include <gsl/gsl_sf_psi.h>
  29. #include <gsl/gsl_sf_gamma.h>
  30.  
  31. #include "error.h"
  32.  
  33. static const double bern[21] = {
  34.    0.0   /* no element 0 */,  
  35.   +0.833333333333333333333333333333333e-01,
  36.   -0.138888888888888888888888888888888e-02,
  37.   +0.330687830687830687830687830687830e-04,
  38.   -0.826719576719576719576719576719576e-06,
  39.   +0.208767569878680989792100903212014e-07,
  40.   -0.528419013868749318484768220217955e-09,
  41.   +0.133825365306846788328269809751291e-10,
  42.   -0.338968029632258286683019539124944e-12,
  43.   +0.858606205627784456413590545042562e-14,
  44.   -0.217486869855806187304151642386591e-15,
  45.   +0.550900282836022951520265260890225e-17,
  46.   -0.139544646858125233407076862640635e-18,
  47.   +0.353470703962946747169322997780379e-20,
  48.   -0.895351742703754685040261131811274e-22,
  49.   +0.226795245233768306031095073886816e-23,
  50.   -0.574472439520264523834847971943400e-24,
  51.   +0.145517247561486490186626486727132e-26,
  52.   -0.368599494066531017818178247990866e-28,
  53.   +0.933673425709504467203255515278562e-30,
  54.   -0.236502241570062993455963519636983e-31
  55. };
  56.  
  57.  
  58. /* ((a)_x - 1)/x in the "small x" region where
  59.  * cancellation must be controlled.
  60.  *
  61.  * Based on SLATEC DPOCH1().
  62.  */
  63. /*
  64. C When ABS(X) is so small that substantial cancellation will occur if
  65. C the straightforward formula is used, we use an expansion due
  66. C to Fields and discussed by Y. L. Luke, The Special Functions and Their
  67. C Approximations, Vol. 1, Academic Press, 1969, page 34.
  68. C
  69. C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as
  70. C        (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .
  71. C In order to maintain significance in POCH1, we write for positive a
  72. C        (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)
  73. C                       = 1.0 + Q*EXPREL(Q) .
  74. C Likewise the polynomial is written
  75. C        POLY = 1.0 + X*POLY1(A,X) .
  76. C Thus,
  77. C        POCH1(A,X) = (POCH(A,X) - 1) / X
  78. C                   = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)
  79. C
  80. */
  81. static
  82. int
  83. pochrel_smallx(const double a, const double x, gsl_sf_result * result)
  84. {
  85.   /*
  86.    SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1))
  87.    ALNEPS = LOG(D1MACH(3))
  88.    */
  89.   const double SQTBIG = 1.0/(2.0*M_SQRT2*M_SQRT3*GSL_SQRT_DBL_MIN);
  90.   const double ALNEPS = GSL_LOG_DBL_EPSILON - M_LN2;
  91.  
  92.   if(x == 0.0) {
  93.     return gsl_sf_psi_e(a, result);
  94.   }
  95.   else {
  96.     const double bp   = (  (a < -0.5) ? 1.0-a-x : a );
  97.     const int     incr = ( (bp < 10.0) ? 11.0-bp : 0 );
  98.     const double b    = bp + incr;
  99.     double dpoch1;
  100.     gsl_sf_result dexprl;
  101.     int stat_dexprl;
  102.     int i;
  103.  
  104.     double var    = b + 0.5*(x-1.0);
  105.     double alnvar = log(var);
  106.     double q = x*alnvar;
  107.  
  108.     double poly1 = 0.0;
  109.  
  110.     if(var < SQTBIG) {
  111.       const int nterms = (int)(-0.5*ALNEPS/alnvar + 1.0);
  112.       const double var2 = (1.0/var)/var;
  113.       const double rho  = 0.5 * (x + 1.0);
  114.       double term = var2;
  115.       double gbern[24];
  116.       int k, j;
  117.  
  118.       gbern[1] = 1.0;
  119.       gbern[2] = -rho/12.0;
  120.       poly1 = gbern[2] * term;
  121.  
  122.       if(nterms > 20) {
  123.         /* NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD */
  124.         /* nterms = 20; */
  125.     result->val = 0.0;
  126.     result->err = 0.0;
  127.     GSL_ERROR ("error", GSL_ESANITY);
  128.       }
  129.  
  130.       for(k=2; k<=nterms; k++) {
  131.         double gbk = 0.0;
  132.         for(j=1; j<=k; j++) {
  133.           gbk += bern[k-j+1]*gbern[j];
  134.         }
  135.         gbern[k+1] = -rho*gbk/k;
  136.  
  137.         term  *= (2*k-2-x)*(2*k-1-x)*var2;
  138.         poly1 += gbern[k+1]*term;
  139.       }
  140.     }
  141.  
  142.     stat_dexprl = gsl_sf_expm1_e(q, &dexprl);
  143.     if(stat_dexprl != GSL_SUCCESS) {
  144.       result->val = 0.0;
  145.       result->err = 0.0;
  146.       return stat_dexprl;
  147.     }
  148.     dexprl.val = dexprl.val/q;
  149.     poly1 *= (x - 1.0);
  150.     dpoch1 = dexprl.val * (alnvar + q * poly1) + poly1;
  151.  
  152.     for(i=incr-1; i >= 0; i--) {
  153.       /*
  154.        C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
  155.        C TO OBTAIN DPOCH1(BP,X).
  156.        */
  157.       double binv = 1.0/(bp+i);
  158.       dpoch1 = (dpoch1 - binv) / (1.0 + x*binv);
  159.     }
  160.  
  161.     if(bp == a) {
  162.       result->val = dpoch1;
  163.       result->err = 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
  164.       return GSL_SUCCESS;
  165.     }
  166.     else {
  167.       /*
  168.        C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A
  169.        C REFLECTION FORMULA TO OBTAIN DPOCH1(A,X).
  170.        */
  171.       double sinpxx = sin(M_PI*x)/x;
  172.       double sinpx2 = sin(0.5*M_PI*x);
  173.       double t1 = sinpxx/tan(M_PI*b);
  174.       double t2 = 2.0*sinpx2*(sinpx2/x);
  175.       double trig  = t1 - t2;
  176.       result->val  = dpoch1 * (1.0 + x*trig) + trig;
  177.       result->err  = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2));
  178.       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(incr) + 1.0) * fabs(result->val);
  179.       return GSL_SUCCESS;
  180.     }    
  181.   }
  182. }
  183.  
  184.  
  185. /* Assumes a>0 and a+x>0.
  186.  */
  187. static
  188. int
  189. lnpoch_pos(const double a, const double x, gsl_sf_result * result)
  190. {
  191.   double absx = fabs(x);
  192.  
  193.   if(absx > 0.1*a || absx*log(GSL_MAX_DBL(a,2.0)) > 0.1) {
  194.     if(a < GSL_SF_GAMMA_XMAX && a+x < GSL_SF_GAMMA_XMAX) {
  195.       /* If we can do it by calculating the gamma functions
  196.        * directly, then that will be more accurate than
  197.        * doing the subtraction of the logs.
  198.        */
  199.       gsl_sf_result g1;
  200.       gsl_sf_result g2;
  201.       gsl_sf_gammainv_e(a,   &g1);
  202.       gsl_sf_gammainv_e(a+x, &g2);
  203.       result->val  = -log(g2.val/g1.val);
  204.       result->err  = g1.err/fabs(g1.val) + g2.err/fabs(g2.val);
  205.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  206.       return GSL_SUCCESS;
  207.     }
  208.     else {
  209.       /* Otherwise we must do the subtraction.
  210.        */
  211.       gsl_sf_result lg1;
  212.       gsl_sf_result lg2;
  213.       int stat_1 = gsl_sf_lngamma_e(a,   &lg1);
  214.       int stat_2 = gsl_sf_lngamma_e(a+x, &lg2);
  215.       result->val  = lg2.val - lg1.val;
  216.       result->err  = lg2.err + lg1.err;
  217.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  218.       return GSL_ERROR_SELECT_2(stat_1, stat_2);
  219.     }
  220.   }
  221.   else if(absx < 0.1*a && a > 15.0) {
  222.     /* Be careful about the implied subtraction.
  223.      * Note that both a+x and and a must be
  224.      * large here since a is not small
  225.      * and x is not relatively large.
  226.      * So we calculate using Stirling for Log[Gamma(z)].
  227.      *
  228.      *   Log[Gamma(a+x)/Gamma(a)] = x(Log[a]-1) + (x+a-1/2)Log[1+x/a]
  229.      *                              + (1/(1+eps)   - 1) / (12 a)
  230.      *                              - (1/(1+eps)^3 - 1) / (360 a^3)
  231.      *                              + (1/(1+eps)^5 - 1) / (1260 a^5)
  232.      *                              - (1/(1+eps)^7 - 1) / (1680 a^7)
  233.      *                              + ...
  234.      */
  235.     const double eps = x/a;
  236.     const double den = 1.0 + eps;
  237.     const double d3 = den*den*den;
  238.     const double d5 = d3*den*den;
  239.     const double d7 = d5*den*den;
  240.     const double c1 = -eps/den;
  241.     const double c3 = -eps*(3.0+eps*(3.0+eps))/d3;
  242.     const double c5 = -eps*(5.0+eps*(10.0+eps*(10.0+eps*(5.0+eps))))/d5;
  243.     const double c7 = -eps*(7.0+eps*(21.0+eps*(35.0+eps*(35.0+eps*(21.0+eps*(7.0+eps))))))/d7;
  244.     const double p8 = gsl_sf_pow_int(1.0+eps,8);
  245.     const double c8 = 1.0/p8             - 1.0;  /* these need not   */
  246.     const double c9 = 1.0/(p8*(1.0+eps)) - 1.0;  /* be very accurate */
  247.     const double a4 = a*a*a*a;
  248.     const double a6 = a4*a*a;
  249.     const double ser_1 = c1 + c3/(30.0*a*a) + c5/(105.0*a4) + c7/(140.0*a6);
  250.     const double ser_2 = c8/(99.0*a6*a*a) - 691.0/360360.0 * c9/(a6*a4);
  251.     const double ser = (ser_1 + ser_2)/ (12.0*a);
  252.  
  253.     double term1 = x * log(a/M_E);
  254.     double term2;
  255.     gsl_sf_result ln_1peps;
  256.     gsl_sf_log_1plusx_e(eps, &ln_1peps);  /* log(1 + x/a) */
  257.     term2 = (x + a - 0.5) * ln_1peps.val;
  258.  
  259.     result->val  = term1 + term2 + ser;
  260.     result->err  = GSL_DBL_EPSILON*fabs(term1);
  261.     result->err += fabs((x + a - 0.5)*ln_1peps.err);
  262.     result->err += fabs(ln_1peps.val) * GSL_DBL_EPSILON * (fabs(x) + fabs(a) + 0.5);
  263.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  264.     return GSL_SUCCESS;
  265.   }
  266.   else {
  267.     gsl_sf_result poch_rel;
  268.     int stat_p = pochrel_smallx(a, x, &poch_rel);
  269.     double eps = x*poch_rel.val;
  270.     int stat_e = gsl_sf_log_1plusx_e(eps, result);
  271.     result->err  = 2.0 * fabs(x * poch_rel.err / (1.0 + eps));
  272.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  273.     return GSL_ERROR_SELECT_2(stat_e, stat_p);
  274.   }
  275. }
  276.  
  277.  
  278. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  279.  
  280. int
  281. gsl_sf_lnpoch_e(const double a, const double x, gsl_sf_result * result)
  282. {
  283.   /* CHECK_POINTER(result) */
  284.  
  285.   if(a <= 0.0 || a+x <= 0.0) {
  286.     DOMAIN_ERROR(result);
  287.   }
  288.   else if(x == 0.0) {
  289.     result->val = 1.0;
  290.     result->err = 0.0;
  291.     return GSL_SUCCESS;
  292.   }
  293.   else {
  294.     return lnpoch_pos(a, x, result);
  295.   }
  296. }
  297.  
  298.  
  299. int
  300. gsl_sf_lnpoch_sgn_e(const double a, const double x,
  301.                        gsl_sf_result * result, double * sgn)
  302. {
  303.   if(a == 0.0 || a+x == 0.0) {
  304.     *sgn = 0.0;
  305.     DOMAIN_ERROR(result);
  306.   }
  307.   else if(x == 0.0) {
  308.     *sgn = 1.0;
  309.     result->val = 1.0;
  310.     result->err = 0.0;
  311.     return GSL_SUCCESS;
  312.   }
  313.   else if(a > 0.0 && a+x > 0.0) {
  314.     *sgn = 1.0;
  315.     return lnpoch_pos(a, x, result);
  316.   }
  317.   else if(a < 0.0 && a+x < 0.0) {
  318.     /* Reduce to positive case using reflection.
  319.      */
  320.     double sin_1 = sin(M_PI * (1.0 - a));
  321.     double sin_2 = sin(M_PI * (1.0 - a - x));
  322.     if(sin_1 == 0.0 || sin_2 == 0.0) {
  323.       *sgn = 0.0;
  324.       DOMAIN_ERROR(result);
  325.     }
  326.     else {
  327.       gsl_sf_result lnp_pos;
  328.       int stat_pp   = lnpoch_pos(1.0-a, -x, &lnp_pos);
  329.       double lnterm = log(fabs(sin_1/sin_2));
  330.       result->val  = lnterm - lnp_pos.val;
  331.       result->err  = lnp_pos.err;
  332.       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(1.0-a) + fabs(1.0-a-x)) * fabs(lnterm);
  333.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  334.       *sgn = GSL_SIGN(sin_1*sin_2);
  335.       return stat_pp;
  336.     }
  337.   }
  338.   else {
  339.     /* Evaluate gamma ratio directly.
  340.      */
  341.     gsl_sf_result lg_apn;
  342.     gsl_sf_result lg_a;
  343.     double s_apn, s_a;
  344.     int stat_apn = gsl_sf_lngamma_sgn_e(a+x, &lg_apn, &s_apn);
  345.     int stat_a   = gsl_sf_lngamma_sgn_e(a,   &lg_a,   &s_a);
  346.     if(stat_apn == GSL_SUCCESS && stat_a == GSL_SUCCESS) {
  347.       result->val  = lg_apn.val - lg_a.val;
  348.       result->err  = lg_apn.err + lg_a.err;
  349.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  350.       *sgn = s_a * s_apn;
  351.       return GSL_SUCCESS;
  352.     }
  353.     else if(stat_apn == GSL_EDOM || stat_a == GSL_EDOM){
  354.       *sgn = 0.0;
  355.       DOMAIN_ERROR(result);
  356.     }
  357.     else {
  358.       result->val = 0.0;
  359.       result->err = 0.0;
  360.       *sgn = 0.0;
  361.       return GSL_FAILURE;
  362.     }
  363.   }
  364. }
  365.  
  366.  
  367. int
  368. gsl_sf_poch_e(const double a, const double x, gsl_sf_result * result)
  369. {
  370.   /* CHECK_POINTER(result) */
  371.  
  372.   if(x == 0.0) {
  373.     result->val = 1.0;
  374.     result->err = 0.0;
  375.     return GSL_SUCCESS;
  376.   }
  377.   else {
  378.     gsl_sf_result lnpoch;
  379.     double sgn;
  380.     int stat_lnpoch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
  381.     int stat_exp    = gsl_sf_exp_err_e(lnpoch.val, lnpoch.err, result);
  382.     result->val *= sgn;
  383.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  384.     return GSL_ERROR_SELECT_2(stat_exp, stat_lnpoch);
  385.   }
  386. }
  387.  
  388.  
  389. int
  390. gsl_sf_pochrel_e(const double a, const double x, gsl_sf_result * result)
  391. {
  392.   const double absx = fabs(x);
  393.   const double absa = fabs(a);
  394.  
  395.   /* CHECK_POINTER(result) */
  396.  
  397.   if(absx > 0.1*absa || absx*log(GSL_MAX(absa,2.0)) > 0.1) {
  398.     gsl_sf_result lnpoch;
  399.     double sgn;
  400.     int stat_poch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);
  401.     if(lnpoch.val > GSL_LOG_DBL_MAX) {
  402.       OVERFLOW_ERROR(result);
  403.     }
  404.     else {
  405.       const double el = exp(lnpoch.val);
  406.       result->val  = (sgn*el - 1.0)/x;
  407.       result->err  = fabs(result->val) * (lnpoch.err + 2.0 * GSL_DBL_EPSILON);
  408.       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(sgn*el) + 1.0) / fabs(x);
  409.       return stat_poch;
  410.     }
  411.   }
  412.   else {
  413.     return pochrel_smallx(a, x, result);
  414.   }
  415. }
  416.  
  417.  
  418. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  419.  
  420. #include "eval.h"
  421.  
  422. double gsl_sf_lnpoch(const double a, const double x)
  423. {
  424.   EVAL_RESULT(gsl_sf_lnpoch_e(a, x, &result));
  425. }
  426.  
  427. double gsl_sf_poch(const double a, const double x)
  428. {
  429.   EVAL_RESULT(gsl_sf_poch_e(a, x, &result));
  430. }
  431.  
  432. double gsl_sf_pochrel(const double a, const double x)
  433. {
  434.   EVAL_RESULT(gsl_sf_pochrel_e(a, x, &result));
  435. }
  436.